function [dof_map_weight, dof_map_weight_inv] = build_c1_weight(V, T, E, TE, d)
%
% function [dof_map_weight, dof_map_weight_inv] = build_c1_weight(V, T, E, TE, d)
% 
% To calculate the d.o.f. weight for the Argyris elements with degree d
%

% calculate the barycentric coordinate of x- and y- axis 
% with respect to all triangles
%
[taux1,taux2,taux3] = directional_bary(V,T,1,0);
[tauy1,tauy2,tauy3] = directional_bary(V,T,0,1);

% calculate the edges' normal direction's barycentric coordinate
nE = size(E,1);
normal = zeros(2,nE);
for eg = 1:nE
    dir = (V(E(eg,2),:) - V(E(eg,1),:))';
    normal(:,eg) = [ 0 -1 ; 1 0]*dir;
%    normal(:,eg) = [ 0 -1 ; 1 0]*(dir/norm(dir,2));
end

% prepare space for dof map weight
m = (d+1)*(d+2)/2;
ndofv = 6; 
ndofe = 2*d - 9;
ndoft = (d-4)*(d-5)/2;
n = 3*ndofv + 3*ndofe + ndoft;

nT = size(T,1);
dof_map_weight = zeros(n,m,nT); 
dof_map_weight_inv = zeros(n,m,nT);

%% 6 d.o.f. per vertex
% vertex 1
    
  dof_map_weight(1,1,:) = 1;
  
  dof_map_weight(2,1,:) = d*taux1; dof_map_weight(3,1,:) = d*tauy1;
  dof_map_weight(2,2,:) = d*taux2; dof_map_weight(3,2,:) = d*tauy2;
  dof_map_weight(2,3,:) = d*taux3; dof_map_weight(3,3,:) = d*tauy3;
  
  dof_map_weight(4,1,:) = d*(d-1)*taux1.*taux1;
  dof_map_weight(4,2,:) = d*(d-1)*2*taux1.*taux2;
  dof_map_weight(4,3,:) = d*(d-1)*2*taux1.*taux3;
  dof_map_weight(4,4,:) = d*(d-1)*taux2.*taux2;
  dof_map_weight(4,5,:) = d*(d-1)*2*taux2.*taux3;
  dof_map_weight(4,6,:) = d*(d-1)*taux3.*taux3;

  dof_map_weight(5,1,:) = d*(d-1)*taux1.*tauy1;
  dof_map_weight(5,2,:) = d*(d-1)*taux1.*tauy2 + d*(d-1)*tauy1.*taux2;
  dof_map_weight(5,3,:) = d*(d-1)*taux1.*tauy3 + d*(d-1)*tauy1.*taux3;
  dof_map_weight(5,4,:) = d*(d-1)*taux2.*tauy2;
  dof_map_weight(5,5,:) = d*(d-1)*taux2.*tauy3 + d*(d-1)*taux3.*tauy2;
  dof_map_weight(5,6,:) = d*(d-1)*taux3.*tauy3;
  
  dof_map_weight(6,1,:) = d*(d-1)*tauy1.*tauy1;
  dof_map_weight(6,2,:) = d*(d-1)*2*tauy1.*tauy2;
  dof_map_weight(6,3,:) = d*(d-1)*2*tauy1.*tauy3;
  dof_map_weight(6,4,:) = d*(d-1)*tauy2.*tauy2;
  dof_map_weight(6,5,:) = d*(d-1)*2*tauy2.*tauy3;
  dof_map_weight(6,6,:) = d*(d-1)*tauy3.*tauy3;
  
% vertex 2

tau = taux1; taux1 = taux2; taux2 = taux3; taux3 = tau;
tau = tauy1; tauy1 = tauy2; tauy2 = tauy3; tauy3 = tau;

  dof_map_weight(7,16,:) = 1;
  
  dof_map_weight(8,16,:) = d*taux1; dof_map_weight(9,16,:) = d*tauy1;
  dof_map_weight(8,17,:) = d*taux2; dof_map_weight(9,17,:) = d*tauy2;
  dof_map_weight(8,11,:) = d*taux3; dof_map_weight(9,11,:) = d*tauy3;
  
  dof_map_weight(10,16,:) = d*(d-1)*taux1.*taux1;
  dof_map_weight(10,17,:) = d*(d-1)*2*taux1.*taux2;
  dof_map_weight(10,11,:) = d*(d-1)*2*taux1.*taux3;
  dof_map_weight(10,18,:) = d*(d-1)*taux2.*taux2;
  dof_map_weight(10,12,:) = d*(d-1)*2*taux2.*taux3;
  dof_map_weight(10,7,:) = d*(d-1)*taux3.*taux3;

  dof_map_weight(11,16,:) = d*(d-1)*taux1.*tauy1;
  dof_map_weight(11,17,:) = d*(d-1)*taux1.*tauy2 + d*(d-1)*tauy1.*taux2;
  dof_map_weight(11,11,:) = d*(d-1)*taux1.*tauy3 + d*(d-1)*tauy1.*taux3;
  dof_map_weight(11,18,:) = d*(d-1)*taux2.*tauy2;
  dof_map_weight(11,12,:) = d*(d-1)*taux2.*tauy3 + d*(d-1)*taux3.*tauy2;
  dof_map_weight(11,7,:) = d*(d-1)*taux3.*tauy3;
  
  dof_map_weight(12,16,:) = d*(d-1)*tauy1.*tauy1;
  dof_map_weight(12,17,:) = d*(d-1)*2*tauy1.*tauy2;
  dof_map_weight(12,11,:) = d*(d-1)*2*tauy1.*tauy3;
  dof_map_weight(12,18,:) = d*(d-1)*tauy2.*tauy2;
  dof_map_weight(12,12,:) = d*(d-1)*2*tauy2.*tauy3;
  dof_map_weight(12,7,:) = d*(d-1)*tauy3.*tauy3;

% vertex 3

tau = taux1; taux1 = taux2; taux2 = taux3; taux3 = tau;
tau = tauy1; tauy1 = tauy2; tauy2 = tauy3; tauy3 = tau;

  dof_map_weight(13,21,:) = 1;
  
  dof_map_weight(14,21,:) = d*taux1; dof_map_weight(15,21,:) = d*tauy1;
  dof_map_weight(14,15,:) = d*taux2; dof_map_weight(15,15,:) = d*tauy2;
  dof_map_weight(14,20,:) = d*taux3; dof_map_weight(15,20,:) = d*tauy3;
  
  dof_map_weight(16,21,:) = d*(d-1)*taux1.*taux1;
  dof_map_weight(16,15,:) = d*(d-1)*2*taux1.*taux2;
  dof_map_weight(16,20,:) = d*(d-1)*2*taux1.*taux3;
  dof_map_weight(16,10,:) = d*(d-1)*taux2.*taux2;
  dof_map_weight(16,14,:) = d*(d-1)*2*taux2.*taux3;
  dof_map_weight(16,19,:) = d*(d-1)*taux3.*taux3;

  dof_map_weight(17,21,:) = d*(d-1)*taux1.*tauy1;
  dof_map_weight(17,15,:) = d*(d-1)*taux1.*tauy2 + d*(d-1)*tauy1.*taux2;
  dof_map_weight(17,20,:) = d*(d-1)*taux1.*tauy3 + d*(d-1)*tauy1.*taux3;
  dof_map_weight(17,10,:) = d*(d-1)*taux2.*tauy2;
  dof_map_weight(17,14,:) = d*(d-1)*taux2.*tauy3 + d*(d-1)*taux3.*tauy2;
  dof_map_weight(17,19,:) = d*(d-1)*taux3.*tauy3;
  
  dof_map_weight(18,21,:) = d*(d-1)*tauy1.*tauy1;
  dof_map_weight(18,15,:) = d*(d-1)*2*tauy1.*tauy2;
  dof_map_weight(18,20,:) = d*(d-1)*2*tauy1.*tauy3;
  dof_map_weight(18,10,:) = d*(d-1)*tauy2.*tauy2;
  dof_map_weight(18,14,:) = d*(d-1)*2*tauy2.*tauy3;
  dof_map_weight(18,19,:) = d*(d-1)*tauy3.*tauy3;

for k = 1:size(T,1)
% %% 2*d - 9 d.o.f. per edge
% % edge 1: v2-v3
%     eg = TE(k,1);
%     [tau1,tau2,tau3] = bary(V(T(k,1),:), V(T(k,2),:), V(T(k,3),:), normal(1,eg), normal(2,eg));
%     dof_map_weight(19,[11 12 13 14 15 16 17 18 19 20 21],k) = [tau1/16 tau1/4 3*tau1/8 tau1/4 tau1/16 ...
%             tau2/16 tau2/4+tau3/16 3*tau2/8+tau3/4 tau2/4+3*tau3/8 tau2/16+tau3/4 tau3/16]; 
% % edge 2: v3-v1
%     eg = TE(k,2);
%     [tau1,tau2,tau3] = bary(V(T(k,1),:), V(T(k,2),:), V(T(k,3),:), normal(1,eg), normal(2,eg));
%     dof_map_weight(20,[20 14 9 5 2 21 15 10 6 3 1],k) = [tau1/16 tau1/4 3*tau1/8 tau1/4 tau1/16 ...
%             tau2/16 tau2/4+tau3/16 3*tau2/8+tau3/4 tau2/4+3*tau3/8 tau2/16+tau3/4 tau3/16]; 
% % edge 3: v1-v2;
%     eg = TE(k,3);
%     [tau1,tau2,tau3] = bary(V(T(k,1),:), V(T(k,2),:), V(T(k,3),:), normal(1,eg), normal(2,eg));
%     dof_map_weight(21,[3 5 8 12 17 1 2 4 7 11 16],k) = [tau1/16 tau1/4 3*tau1/8 tau1/4 tau1/16 ...
%             tau2/16 tau2/4+tau3/16 3*tau2/8+tau3/4 tau2/4+3*tau3/8 tau2/16+tau3/4 tau3/16]; 

%% There is no identity matrix for the rest inner d.o.f.

    idx = [1 2 3 4 5 6 16 17 11 18 12 7 21 15 20 10 14 19 8 9 13];
    dof_map_weight_inv(idx,:,k) = inv(dof_map_weight(:,idx,k));

end

end

function [tau1,tau2,tau3] = bary(V1,V2,V3,x,y)
x1 = V1(1); x2 = V2(1); x3 = V3(1);
y1 = V1(2); y2 = V2(2); y3 = V3(2);

J  = (x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1);
J1 = (x2 - x )*(y3 - y ) - (x3 - x )*(y2 - y );
J2 = (x  - x1)*(y3 - y1) - (x3 - x1)*(y  - y1);
J3 = (x2 - x1)*(y  - y1) - (x  - x1)*(y2 - y1);

% barycetric for (0,0);
J10 = x2.*y3 - x3.*y2;
J20 = (x3 - x1).*y1 - (y3 - y1).*x1;
J30 = (y2 - y1).*x1 - (x2 - x1).*y1;

tau1 = (J1 - J10)./J;
tau2 = (J2 - J20)./J;
tau3 = (J3 - J30)./J;
end


